
####################以下内容需要根据自己的需要进行修改#####################

BigAlgorithm=function(outpath,gen,myY,method,myYs){
	#设置工作目录
	#outpath="C:\\Users\\Mr-LiZhuo\\Desktop\\13种算法\\11"
	#setwd(outpath)
	#读取文件，其中二者必须有个体名称
	#gen=read.table("../mdp_numeric.txt",head=F)#基因型文件，基因型文件的第一列为个体id，从第二列开始往后是基因值
	#myY=read.table("../mdp_Y.txt",head=F)#表型文件，表型文件的第一列为个体id，第二列为表型值
	#设定计算的方法，方便提取结果
	#method=c("RandomForest","SVM","BayesA","BayesB","BayesC","BayesL","BRR","RKHS","LASSO","RidgeReg","ElasticNet","PLS","BayesG","GBLUP")
	#myYs=c(2:4,5)
	setwd(outpath)	

	############################以下内容可以不修改############################
	#导入必须的计算包
	if(!require("BGLR")){
		install.packages("BGLR")
	}else{
		require("BGLR",quietly = T)
	}

	if(!require("ranger")){
		install.packages("ranger")
	}else{
		require("ranger",quietly = T)
	}

	if(!require("kernlab")){
		install.packages("kernlab")
	}else{
		require("kernlab",quietly = T)
	} 
	
	if(!require("glmnet")){
		install.packages("glmnet")
	}else{
		require("glmnet",quietly = T)
	}

	if(!require("NAM")){
		install.packages("NAM")
	}else{
		require("NAM",quietly = T)
	}

	if(!require("rrBLUP")){
		install.packages("rrBLUP")
	}else{
		require("rrBLUP",quietly = T)
	}

	if(!require("pls")){
		install.packages("pls")
	}else{
		require("pls",quietly = T)
	}
	
	#将基因型文件和表型文件的顺讯进行匹配，匹配的结果以表型文件为准
	myY<-myY[myY[,1] %in% gen[,1],]
	gen<-gen[gen[,1] %in% myY[,1],]
	order=match(myY[,1],gen[,1])
	gen=gen[order,]	

	#去除基因型文件的行标题，方便后面计算
	gen = gen[,-1]	

	#将gen文件转化为matrix格式的，方便后面计算
	gen=as.matrix(gen)	

	###################################注意，随机森林算法非常费内存，运行速度非常慢，如果设备性能不足，请屏蔽随机森林算法########################	

	for(nrows in 1:length(myYs)){#按列循环
		y=as.numeric(unlist(myY[,myYs[nrows]]))#将表型文件转化为numeric格式，方便后面计算
		w=which(is.na(y))#提取出NA所在的行坐标
		for(i in 1:length(method)){#按方法循环
			if(method[i]=="RandomForest"){
				#RandomForest随机森林算法		
				if(length(w)!=0){
					temp1 = ranger(y~.,data.frame(y=y[-w],gen[-w,]),importance='impurity',num.trees=1000)
				}else{
					temp1 = ranger(y~.,data.frame(y,gen),importance='impurity',num.trees=1000)
				}		
				RandomForest = predict(temp1,data.frame(gen))$predictions
				write.table(cbind(myY[,1],RandomForest),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    
			}else if(method[i]=="SVM"){
				#SVM算法
				SVM = predict(ksvm(gen,y,type= "eps-svr"),gen)[,1]; 
				write.table(cbind(myY[,1],SVM),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="BayesA"){
				#Bayesian系列算法
				BayesA = BGLR(rmExistingFiles = F,y,ETA=list(list(X=gen,model='BayesA')),verbose=F)$yHat
				write.table(cbind(myY[,1],BayesA),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="BayesB"){
				#Bayesian系列算法
				BayesB = BGLR(rmExistingFiles = F,y,ETA=list(list(X=gen,model='BayesB')),verbose=F)$yHat
				write.table(cbind(myY[,1],BayesB),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="BayesC"){
				#Bayesian系列算法
				BayesC = BGLR(rmExistingFiles = F,y,ETA=list(list(X=gen,model='BayesC')),verbose=F)$yHat
				write.table(cbind(myY[,1],BayesC),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="BayesL"){
				#Bayesian系列算法
				BayesL = BGLR(rmExistingFiles = F,y,ETA=list(list(X=gen,model='BL')),verbose=F)$yHat
				write.table(cbind(myY[,1],BayesL),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="BRR"){
				#Bayesian系列算法
				BRR = BGLR(rmExistingFiles = F,y,ETA=list(list(X=gen,model='BRR')),verbose=F)$yHat
				write.table(cbind(myY[,1],BRR),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="RKHS"){
				#Bayesian系列算法
				#计算G矩阵，方便GBLUP、BayesG、RKHS三个方法计算。如果没有用这三种方法，那么以下四行代码可以不运行。
				E2 = as.matrix(dist(gen)^2)
				G = exp(-(E2/mean(E2)))
				eK = eigen(G,symmetric=T)
				RKHS = BGLR(rmExistingFiles = F,y,ETA=list(list(V=eK$vectors,d=eK$values,model='RKHS')),verbose=F)$yHat
				write.table(cbind(myY[,1],RKHS),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="LASSO"){
				#LASSO算法	
				if(length(w)!=0){
					cv1 = glmnet::cv.glmnet(x=gen[-w,],y=y[-w],alpha=1); 
					lmb1 = cv1$lambda.min;
					f4 = glmnet::glmnet(x=gen[-w,],y=y[-w],lambda = lmb1,alpha = 1)
				}else{
					cv1 = glmnet::cv.glmnet(x=gen,y=y,alpha=1); 
					lmb1 = cv1$lambda.min;
					f4 = glmnet::glmnet(x=gen,y=y,lambda = lmb1,alpha = 1)
				}
				LASSO = c(glmnet::predict.glmnet(f4,gen))
				write.table(cbind(myY[,1],LASSO),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="RidgeReg"){
				#RidgeReg算法	
				if(length(w)!=0){
					cv2 = glmnet::cv.glmnet(x=gen[-w,],y=y[-w],alpha=0);
					lmb2 = cv2$lambda.min; 
					f5 = glmnet::glmnet(x=gen[-w,],y=y[-w],lambda = lmb2,alpha = 0)		
				}else{
					cv2 = glmnet::cv.glmnet(x=gen,y,alpha=0);
					lmb2 = cv2$lambda.min; 
					f5 = glmnet::glmnet(x=gen,y=y,lambda = lmb2,alpha = 0)
				}
				RidgeReg = c(glmnet::predict.glmnet(f5,gen))
				write.table(cbind(myY[,1],RidgeReg),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="ElasticNet"){
				#ElasticNet算法	
				if(length(w)!=0){
					cv3 = glmnet::cv.glmnet(x=gen[-w,],y=y[-w],alpha=0.01);
					lmb3 = cv3$lambda.min; 
					f14 = glmnet::glmnet(x=gen[-w,],y=y[-w],lambda = lmb3,alpha = 0.01)		
				}else{
					cv3 = glmnet::cv.glmnet(x=gen,y,alpha=0.01);
					lmb3 = cv3$lambda.min; 
					f14 = glmnet::glmnet(x=gen,y,lambda = lmb3,alpha = 0.01)
				}
				ElasticNet = c(glmnet::predict.glmnet(f14,gen))
				write.table(cbind(myY[,1],ElasticNet),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="PLS"){
				#PLS算法	
				if(length(w)!=0){
					ncomp = 5
					f16a = plsr(y[-w]~gen[-w,],ncomp=ncomp,validation='none')				
				}else{
					ncomp = 5
					f16a = plsr(y~gen,ncomp=ncomp,validation='none')		
				}	
				PLS = predict(f16a,gen)[,1,ncomp]	
				0
			}else if(method[i]=="BayesG"){
				#BayesG算法（RKHS算法采用G矩阵计算）
				G = NAM::GRM(gen,TRUE)
				BayesG=BGLR( y=y,ETA=list(G=list(K=G,model='RKHS')))$yHat	
				write.table(cbind(myY[,1],BayesG),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}else if(method[i]=="GBLUP"){
				#GBLUP算法
				G = NAM::GRM(gen,TRUE)	
				f10 = mixed.solve(y,K=G)	
				GBLUP= c(f10$u)	
				write.table(cbind(myY[,1],GBLUP),paste(myYs[nrows],"_",method[i],"_result.txt",sep=''),quote=F,row.names=F,col.names = F)    	

			}		
		}	
	}
}
